Part 2 Fire Selection
2.1 Set Up
2.1.1 Libraries
library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(patchwoRk)
library(gridExtra)2.1.2 USDA National Forest Type Group Dataset
Conifer Forest Type Groups: Douglas-Fir, Fir-Spruce-Mountain Hemlock, Lodgepole Pine
# forest type groups and key
conus_forestgroup <- raster('data/forest_type/conus_forestgroup.tif')
forest_codes <- read_csv('data/forest_type/forestgroupcodes.csv')
# set crs
crs = crs(conus_forestgroup)2.1.3 EPA level-3 Ecoregions
Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies
# level 3 ecoregions
l3eco <- st_read('data/ecoregion/us_eco_l3.shp') %>%
st_transform(., crs=crs)
# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>%
filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))2.2 Define Fire Parameters
2.2.1 Monitoring Trends in Burn Severity (MTBS) Dataset
Criteria:
-1988-1991
-500+ acres of high-severity
-Within selected ecoregions
->25% of selected forest types
# mtbs fire perimeters
mtbs_full <- st_read('data/mtbs/mtbs_perims_DD.shp') %>%
st_transform(., crs=crs)
mtbs_select <- mtbs_full %>%
mutate(state = str_sub(Event_ID,0,2),
year = year(as.Date(Ig_Date))) %>%
filter(state %in% c("WA","ID","MT","WY","SD"),
between(Ig_Date, as.Date('1988-01-1'), as.Date('1991-12-31'))) 2.2.2 Group Adjacent Fires
# function to group adjoining fire polygons to ensure contiguous high-severity patches
group_fires <- function(mtbs_year) {
# join the polygons with themselves, and remove those that do not join with any besides themselves
combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>%
drop_na(Event_ID.y)%>%
dplyr::select(Event_ID.x,Event_ID.y)
if(nrow(combined)>=1){ # if there are overlaps for this years fires...
# partition data into that that has overlap, and that that does not
overlap <- mtbs_year %>%
filter(Event_ID %in% combined$Event_ID.x)
no_overlap <- mtbs_year %>%
filter(!(Event_ID %in% combined$Event_ID.x))
print(paste0("there are ",nrow(overlap)," overlapping polygons"))
# join all overlapping features, and buffer to ensure proper grouping
overlap_union <- st_union(overlap) %>%
st_buffer(190)
# break apart the joined polygons into their individual groups
groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
mutate(year = mean(mtbs_year$year),
Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
rename(geometry = x)
print(paste0("polygons formed into ",nrow(groups)," groups"))
# join back with original dataset to return to unbuffered geometry
grouped_overlap <- st_join(overlap,groups,left=TRUE)
# arrange by the new grouping
joined_overlap_groups <- grouped_overlap %>%
group_by(Fire_ID) %>%
tally()%>%
st_buffer(1) %>%
dplyr::select(Fire_ID) %>%
mutate(year = mean(mtbs_year$year))
# add new ID to the freestanding polygons
no_overlap_groups <- no_overlap %>%
mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
dplyr::select(Fire_ID,year)
# join the new grouped overlap and the polygons without overlap
fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
return(fires_export)
} else { # if there are no overlaps for this year...
print("no overlapping polygons")
fires_export <- mtbs_year %>%
mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
dplyr::select(Fire_ID,year)
return(fires_export)
}
}# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>% filter(year == 1988))## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"
fires_89 <- group_fires(mtbs_select %>% filter(year == 1989))## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_90 <- group_fires(mtbs_select %>% filter(year == 1990))## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_91 <- group_fires(mtbs_select %>% filter(year == 1991))## [1] "no overlapping polygons"
# join each fire year, filter by area
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
mutate(area_ha = as.numeric(st_area(geometry))/10000,
area_acres = area_ha*2.471)2.2.3 Select Fires by Ecoregion and Forest Type
# assign ecoregion and proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>%
left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08,
fun = function(value, coverage_fraction) {
data.frame(value = value,
frac = coverage_fraction / sum(coverage_fraction)) %>%
group_by(value) %>%
summarize(freq = sum(frac), .groups = 'drop') %>%
pivot_wider(names_from = 'value',
names_prefix = 'freq_',
values_from = 'freq')}) %>%
mutate(across(starts_with('freq'), replace_na, 0)))
# remove unnecessary columns, cleanup names
# filter to ensure fire polygons are at least 25% type of interest
fires <- fires_join %>%
dplyr::select("Fire_ID","year","area_ha","area_acres","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>%
rename("ecoregion" = "US_L3NAME",
"freq_df"="freq_200",
"freq_pp"="freq_220",
"freq_fs"="freq_260",
"freq_lpp"="freq_280") %>%
mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
freq_forested = 1- freq_0,
freq_ideal = freq_df+freq_fs+freq_lpp)%>%
mutate(across(starts_with('freq'), round,2))%>%
filter(freq_ideal > 0.25)2.2.4 Select Fires by Burn Severity
# import all mtbs rasters via a list
rastlist <- list.files(path = "data/mtbs", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,22,25))
# create empty dataframe
severity_list <- list()
# loop through mtbs mosasics for 1988-1991
# extract mtbs burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
mtbs_year <- allrasters[[i]]
fire_year <- filter(fires, year==str_sub(i,2,5))
raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
names(raster_extract) <- fire_year$Fire_ID
output_select <- bind_rows(raster_extract, .id = "Fire_ID")%>%
group_by(Fire_ID , value) %>%
summarize(total_area = sum(coverage_area)) %>%
group_by(Fire_ID) %>%
mutate(proportion = total_area/sum(total_area))%>%
dplyr::select("Fire_ID","value","proportion") %>%
spread(.,key="value",value = "proportion")
severity_list[[i]] <- output_select
}
# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list)
# join burn severity % to fires polygons
# fix naming
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="Fire_ID")%>%
rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>%
dplyr::select(- "NaN",-"regrowth",-"error") %>%
mutate(highsev_acres = area_acres*highsev)%>%
filter(highsev_acres > 500)2.2.5 Clean Up Dataset
# get the most common forest type within each polygon
fires_select <- fires_severity %>%
left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08))
fires_select$mode <- as.factor(fires_select$mode)
fires_select <- fires_select %>%
mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
mode==220 ~ "Ponderosa",
mode==260 ~ "Fir-Spruce",
mode==280 ~ "Lodegepole Pine",
TRUE ~ "Other"),
Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year))# join the grouped fires back to original mtbs boundaries
fires_mtbs <- st_join(mtbs_select,fires_select,left=FALSE,largest=TRUE) %>%
filter(year.x==year.y)%>%
dplyr::select("Event_ID","Incid_Name","Fire_ID","Ig_Date","year.y","state","BurnBndAc","ecoregion") %>%
rename(year= year.y)